###############################################################
### Replication Files for:                                  ###
### Western Political Rhetoric and Radicalization           ###
### William O'Brochta, Margit Tavits, and Deniz Aksoy       ###
### Washington University in St. Louis                      ###
### Published in the British Journal of Political Science   ###
###############################################################

library(foreign)
library(readstata13)
library(haven)
library(AER)
library(interplot)
library(arm)
library(stargazer)
library(xtable)
library(nnet)
library(coefplot2)


#Preparing Dataset
data<-read_dta('MergedDataset.dta')
data$tr_IslamPositive<-ifelse(is.na(data$g2),0,1)
data$tr_IslamNegative<-ifelse(is.na(data$g3),0,1)
data$dv_USfavorable<-ifelse(data$g4==1,1,0)
data$dv_Muslimidentity<-6-data$g5
data$dv_UnderstandViolence<-6-data$g6

#Recode NAs
data$unempl<-ifelse(data$dm7==888997, NA, data$dm7)
data$dm2<-ifelse(data$dm2==888999, NA, data$dm2)
data$dm7<-ifelse(data$dm7==888997, NA, data$dm7)
data$dm24<-ifelse(data$dm24==888999, NA, data$dm24)
data$dm34<-ifelse(data$dm34==888998, NA,data$dm34)

#Month variable
data$month<-ifelse(data$april==1, 2, NA)
data$month<-ifelse(data$may==1, 3, data$month)
data$month<-ifelse(data$june==1, 4, data$month)
data$month<-ifelse(data$july==1, 5, data$month)
data$month<-ifelse(is.na(data$month), 1, data$month)
data$month<-as.factor(data$month)

#unemployed
data$unempl<-ifelse(!is.na(data$unempl) & data$unempl==5, 1,0)
#edu includes university education only
data$edu<-ifelse(data$dm3_edu==3, 1,0)
data$dv_UnderstandViolence<-as.factor(data$dv_UnderstandViolence)
data$dv_Muslimidentity<-as.factor(data$dv_Muslimidentity)
data$edu_factor<-as.factor(data$dm3_edu)
data$dv_USfavorable<-as.factor(data$dv_USfavorable)
data$married<-ifelse(data$dm19==1, 1,0)

#Create control variable
data$control<-ifelse(data$tr_IslamNegative==0 & data$tr_IslamPositive==0, 1,0)
data$assignment<-ifelse(data$control==1, 0, NA)
data$assignment<-ifelse(data$tr_IslamPositive==1, 1, data$assignment)
data$assignment<-ifelse(data$tr_IslamNegative==1, 2, data$assignment)

#SI.5: Comparison of Means
#Muslim Identity
t.test(as.numeric(data[data$tr_IslamNegative==1,]$dv_Muslimidentity),
       as.numeric(data[data$control==1,]$dv_Muslimidentity))
t.test(as.numeric(data[data$tr_IslamPositive==1,]$dv_Muslimidentity),
       as.numeric(data[data$control==1,]$dv_Muslimidentity))

mean(as.numeric(data[data$tr_IslamNegative==1,]$dv_Muslimidentity))
mean(as.numeric(data[data$control==1,]$dv_Muslimidentity))
mean(as.numeric(data[data$tr_IslamPositive==1,]$dv_Muslimidentity))

#U.S. Favorable (coded here as 1=unfavorable, 2=favorable; 
#discussed in SI as percent favorable i.e., subtract 1 from these results)
t.test(as.numeric(data[data$tr_IslamNegative==1,]$dv_USfavorable),
       as.numeric(data[data$control==1,]$dv_USfavorable))
t.test(as.numeric(data[data$tr_IslamPositive==1,]$dv_USfavorable),
       as.numeric(data[data$control==1,]$dv_USfavorable))

mean(as.numeric(data[data$tr_IslamNegative==1,]$dv_USfavorable))
mean(as.numeric(data[data$control==1,]$dv_USfavorable))
mean(as.numeric(data[data$tr_IslamPositive==1,]$dv_USfavorable))

#Understand Violence
t.test(as.numeric(data[data$tr_IslamNegative==1,]$dv_UnderstandViolence),
       as.numeric(data[data$control==1,]$dv_UnderstandViolence))
t.test(as.numeric(data[data$tr_IslamPositive==1,]$dv_UnderstandViolence),
       as.numeric(data[data$control==1,]$dv_UnderstandViolence))

mean(as.numeric(data[data$tr_IslamNegative==1,]$dv_UnderstandViolence))
mean(as.numeric(data[data$control==1,]$dv_UnderstandViolence))
mean(as.numeric(data[data$tr_IslamPositive==1,]$dv_UnderstandViolence))




#Main Text Figure 1 panels a, b, c
source('CoefPlot.R')
pdf("dv_MuslimIdentity.pdf", width=8.5, height=11)
par(mar = c(0, 12, 0, 0))
#Figure 1: Muslim Identity
#dv_MuslimIdentity
coefplot2(c(0.06238,-0.04859),c(0.06063,0.06121), 
         varnames=c('Pro-Islam', 'Anti-Islam'), xlim=c(-.2,.4), main='', 
         mar=c(12,9,3,1), cex.var=2.5, cex.pts=2, h.axis=F)
axis(3,cex.axis=1.7)
dev.off()

pdf("dv_USfavorable.pdf", width=8.5, height=11)
#Figure 2: US Favorable
#dv_USfavorable
coefplot2(c(0.061467,0.006239),c(0.023685,0.023912), 
          varnames=c('', ''), xlim=c(-.2,.4), main='', 
          mar=c(12,9,3,1), cex.var=2.5, cex.pts=2, h.axis=F)
axis(3,cex.axis=1.7)
dev.off()

pdf("dv_UnderstandViolence.pdf", width=8.5, height=11)
#Figure 3: Approve Violence
#dv_UnderstandViolence
coefplot2(c(0.09761,0.07425),c(0.06179,0.06238), 
         varnames=c('', ''), xlim=c(-.2,.4), main='', 
         mar=c(12,9,3,1), cex.var=2.5, cex.pts=2, h.axis=F)
axis(3,cex.axis=1.7)
dev.off()




#SI.2: Additional Information about the Survey
data$immigration<-rowSums(data[,c(3:5)], na.rm=T)
data$immigration<-as.factor(data$immigration)
model<-polr(immigration~tr_IslamPositive+
                 tr_IslamNegative+dm1+dm2+dm3_edu+unempl+married
               +april+may+june+july, 
               data=data, method=c('probit'))

#Table SI.2.1: Immigration Prime
stargazer(model, star.cutoffs = c(0.05, 0.01, 0.001))


#SI.3: Randomization and Balance Checks
#Gender, age, education, employment status, marrital status
#Individual covariates

#Randomization Check
model0<-multinom(assignment~dm1+dm2+dm3_edu+unempl+married, data=data)
z <- summary(model0)$coefficients/summary(model0)$standard.errors
p <- (1 - pnorm(abs(z), 0, 1)) * 2
beta <- as.vector(t(coef(model0)))
A<-rbind(c(1,0,0,0,0,0,-1,0,0,0,0,0), c(0,1,0,0,0,0,0,-1,0,0,0,0),
         c(0,0,1,0,0,0,0,0,-1,0,0,0), c(0,0,0,1,0,0,0,0,0,-1,0,0),
         c(0,0,0,0,1,0,0,0,0,0,-1,0), c(0,0,0,0,0,1,0,0,0,0,0,-1))
t(A %*% beta) %*% solve(A %*% vcov(model0) %*% t(A), A %*% beta) 
pchisq(4.130906, 5, lower=F)

#Table SI.3.1: Randomization Check
stargazer(model0, star.cutoffs = c(0.05, 0.01, 0.001))


#Table SI.3.2: Balance Checks
#Tr 1: Islam Positive
model0.1<-lm(tr_IslamPositive~dm1, data=data)
mean(data[data$tr_IslamPositive==1,]$dm1)
mean(data[data$tr_IslamPositive==0,]$dm1)
model0.2<-lm(tr_IslamPositive~dm2, data=data)
mean(data[data$tr_IslamPositive==1,]$dm2)
mean(data[data$tr_IslamPositive==0,]$dm2)
model0.3<-lm(tr_IslamPositive~dm3_edu, data=data)
mean(data[data$tr_IslamPositive==1,]$dm3_edu)
mean(data[data$tr_IslamPositive==0,]$dm3_edu)
model0.4<-lm(tr_IslamPositive~unempl, data=data)
mean(data[data$tr_IslamPositive==1,]$unempl, na.rm=T)
mean(data[data$tr_IslamPositive==0,]$unempl, na.rm=T)
model0.5<-lm(tr_IslamPositive~married, data=data)
mean(data[data$tr_IslamPositive==1,]$married)
mean(data[data$tr_IslamPositive==0,]$married)

#Tr 2: Islam Negative
model0.7<-lm(tr_IslamNegative~dm1, data=data)
mean(data[data$tr_IslamNegative==1,]$dm1)
mean(data[data$tr_IslamNegative==0,]$dm1)
model0.8<-lm(tr_IslamNegative~dm2, data=data)
mean(data[data$tr_IslamNegative==1,]$dm2)
mean(data[data$tr_IslamNegative==0,]$dm2)
model0.9<-lm(tr_IslamNegative~dm3_edu, data=data)
mean(data[data$tr_IslamNegative==1,]$dm3_edu)
mean(data[data$tr_IslamNegative==0,]$dm3_edu)
model0.10<-lm(tr_IslamNegative~unempl, data=data)
mean(data[data$tr_IslamNegative==1,]$unempl, na.rm=T)
mean(data[data$tr_IslamNegative==0,]$unempl, na.rm=T)
model0.11<-lm(tr_IslamNegative~married, data=data)
mean(data[data$tr_IslamNegative==1,]$married)
mean(data[data$tr_IslamNegative==0,]$married)

#Control
model0.13<-lm(control~dm1, data=data)
mean(data[data$control==1,]$dm1)
mean(data[data$control==0,]$dm1)
model0.14<-lm(control~dm2, data=data)
mean(data[data$control==1,]$dm2)
mean(data[data$control==0,]$dm2)
model0.15<-lm(control~dm3_edu, data=data)
mean(data[data$control==1,]$dm3_edu)
mean(data[data$control==0,]$dm3_edu)
model0.16<-lm(control~unempl, data=data)
mean(data[data$control==1,]$unempl, na.rm=T)
mean(data[data$control==0,]$unempl, na.rm=T)
model0.17<-lm(control~married, data=data)
mean(data[data$control==1,]$married)
mean(data[data$control==0,]$married)

#Wald Test
model0.18<-lm(tr_IslamPositive~dm1+dm2+dm3_edu+unempl+married, data=data)
waldtest(model0.18)

model0.19<-lm(tr_IslamNegative~dm1+dm2+dm3_edu+unempl+married, data=data)
waldtest(model0.19)

model0.20<-lm(control~dm1+dm2+dm3_edu+unempl+married, data=data)
waldtest(model0.20)


#SI.4: Regression Output
model1<-polr(dv_UnderstandViolence~tr_IslamPositive+tr_IslamNegative+
               april+may+june+july, 
             data=data, method=c('probit'))

data$dv_UnderstandViolenceN<-as.numeric(as.character(data$dv_UnderstandViolence))
model1a<-lm(dv_UnderstandViolenceN~tr_IslamPositive+tr_IslamNegative+
              april+may+june+july, data=data)

model1.1<-polr(dv_Muslimidentity~tr_IslamPositive+tr_IslamNegative+
                 april+may+june+july, 
               data=data, method=c('probit'))

data$dv_MuslimidentityN<-as.numeric(as.character(data$dv_Muslimidentity))
model1.1a<-lm(dv_MuslimidentityN~tr_IslamPositive+tr_IslamNegative+
              april+may+june+july, data=data)

model1.2<-glm(dv_USfavorable~tr_IslamPositive+tr_IslamNegative+
               april+may+june+july, 
             data=data, family = binomial(link = "probit"))

data$dv_USfavorableN<-as.numeric(as.character(data$dv_USfavorable))
model1.2a<-lm(dv_USfavorableN~tr_IslamPositive+tr_IslamNegative+
                april+may+june+july, data=data)

#Table SI.4.1: Full Sample With Month Fixed Effects
stargazer(model1.1, model1.1a, model1.2, model1.2a, model1, model1a, star.cutoffs = c(0.05, 0.01, 0.001))


#SI.5: Robustness Checks
model2.1<-polr(dv_UnderstandViolence~tr_IslamPositive+
               tr_IslamNegative+dm1+dm2+dm3_edu+unempl+married
               +april+may+june+july, 
             data=data, method=c('probit'))

model2.2<-polr(dv_Muslimidentity~tr_IslamPositive+tr_IslamNegative+
                 dm1+dm2+dm3_edu+unempl+married
               +april+may+june+july, 
               data=data, method=c('probit'))

model2.3<-glm(dv_USfavorable~tr_IslamPositive+tr_IslamNegative+
                dm1+dm2+dm3_edu+unempl+married
              +april+may+june+july, 
              data=data, family = binomial(link = "probit"))

#March Only
model2.4<-polr(dv_UnderstandViolence~tr_IslamPositive+
               tr_IslamNegative+dm1+dm2+dm3_edu+unempl+married, 
             data=data[data$month==1,], method=c('probit'))
model2.5<-polr(dv_Muslimidentity~tr_IslamPositive+tr_IslamNegative+
                 dm1+dm2+dm3_edu+unempl+married, 
               data=data[data$month==1,], method=c('probit'))
model2.6<-glm(dv_USfavorable~tr_IslamPositive+tr_IslamNegative+
                dm1+dm2+dm3_edu+unempl+married, 
              data=data[data$month==1,], family = binomial(link = "probit"))

#April Only
model2.7<-polr(dv_UnderstandViolence~tr_IslamPositive+
               tr_IslamNegative+dm1+dm2+dm3_edu+unempl+married, 
             data=data[data$april==1,], method=c('probit'))
model2.8<-polr(dv_Muslimidentity~tr_IslamPositive+tr_IslamNegative+
                 dm1+dm2+dm3_edu+unempl+married, 
               data=data[data$april==1,], method=c('probit'))
model2.9<-glm(dv_USfavorable~tr_IslamPositive+tr_IslamNegative+
                dm1+dm2+dm3_edu+unempl+married, 
              data=data[data$april==1,], family = binomial(link = "probit"))

#May Only
model2.10<-polr(dv_UnderstandViolence~tr_IslamPositive+
               tr_IslamNegative+dm1+dm2+dm3_edu+unempl+married, 
             data=data[data$may==1,], method=c('probit'))
model2.11<-polr(dv_Muslimidentity~tr_IslamPositive+tr_IslamNegative+
                 dm1+dm2+dm3_edu+unempl+married, 
               data=data[data$may==1,], method=c('probit'))
model2.12<-glm(dv_USfavorable~tr_IslamPositive+tr_IslamNegative+
                dm1+dm2+dm3_edu+unempl+married, 
              data=data[data$may==1,], family = binomial(link = "probit"))

#June Only
model2.13<-polr(dv_UnderstandViolence~tr_IslamPositive+
               tr_IslamNegative+dm1+dm2+dm3_edu+unempl+married, 
             data=data[data$june==1,], method=c('probit'))
model2.14<-polr(dv_Muslimidentity~tr_IslamPositive+tr_IslamNegative+
                 dm1+dm2+dm3_edu+unempl+married, 
               data=data[data$june==1,], method=c('probit'))
model2.15<-glm(dv_USfavorable~tr_IslamPositive+tr_IslamNegative+
                dm1+dm2+dm3_edu+unempl+married, 
              data=data[data$june==1,], family = binomial(link = "probit"))

#July Only
model2.16<-polr(dv_UnderstandViolence~tr_IslamPositive+
               tr_IslamNegative+dm1+dm2+dm3_edu+unempl+married, 
             data=data[data$july==1,], method=c('probit'))
model2.17<-polr(dv_Muslimidentity~tr_IslamPositive+tr_IslamNegative+
                 dm1+dm2+dm3_edu+unempl+married, 
               data=data[data$july==1,], method=c('probit'))
model2.18<-glm(dv_USfavorable~tr_IslamPositive+tr_IslamNegative+
                dm1+dm2+dm3_edu+unempl+married, 
              data=data[data$july==1,], family = binomial(link = "probit"))

#Table SI.5.1: Muslim Identity Full Sample By Month
stargazer(model2.2, model2.5, model2.8, 
          model2.11, model2.14, model2.17,
          star.cutoffs = c(0.05, 0.01, 0.001))

#Table SI.5.2: U.S. Favorable Full Sample By Month
stargazer(model2.3, model2.6, model2.9,
          model2.12, model2.15, model2.18, 
          star.cutoffs = c(0.05, 0.01, 0.001))

#Table SI.5.3: Approve Violence Full Sample By Month
stargazer(model2.1, model2.4, model2.7,
          model2.10, model2.13, model2.16, 
          star.cutoffs = c(0.05, 0.01, 0.001))



#SI.6: Observational Analysis
#Hegerogeneous Treatment Effects
#Gender
model4.1<-polr(dv_UnderstandViolence~tr_IslamPositive*dm1+
                 tr_IslamNegative*dm1+dm1+dm2+dm3_edu+unempl+married
               +april+may+june+july, 
               data=data, method=c('probit'))
coeftest(model4.1)

model4.2<-polr(dv_Muslimidentity~tr_IslamPositive*dm1+tr_IslamNegative*dm1+
                 dm1+dm2+dm3_edu+unempl+married
               +april+may+june+july, 
               data=data, method=c('probit'))
coeftest(model4.2)

model4.13<-glm(dv_USfavorable~tr_IslamPositive*dm1+tr_IslamNegative*dm1+
                 dm1+dm2+dm3_edu+unempl+married
               +april+may+june+july, 
               data=data, family = binomial(link = "probit"))
coeftest(model4.13)


#Age
model4.3<-polr(dv_UnderstandViolence~tr_IslamPositive*dm2+
                 tr_IslamNegative*dm2+dm1+dm2+dm3_edu+unempl+married
               +april+may+june+july, 
               data=data, method=c('probit'))
coeftest(model4.3)

model4.4<-polr(dv_Muslimidentity~tr_IslamPositive*dm2+tr_IslamNegative*dm2+
                 dm1+dm2+dm3_edu+unempl+married
               +april+may+june+july, 
               data=data, method=c('probit'))
coeftest(model4.4)

model4.14<-glm(dv_USfavorable~tr_IslamPositive*dm2+tr_IslamNegative*dm2+
                 dm1+dm2+dm3_edu+unempl+married
               +april+may+june+july, 
               data=data, family = binomial(link = "probit"))
coeftest(model4.14)


#Married
model4.7<-polr(dv_UnderstandViolence~tr_IslamPositive*married+
                 tr_IslamNegative*married+dm3_edu+unempl+dm1+dm2
               +april+may+june+july, 
               data=data, method=c('probit'))
coeftest(model4.7)

model4.8<-polr(dv_Muslimidentity~tr_IslamPositive*married+
                 tr_IslamNegative*married+dm3_edu+unempl+dm1+dm2
               +april+may+june+july, 
               data=data, method=c('probit'))
coeftest(model4.8)

model4.16<-glm(dv_USfavorable~tr_IslamPositive*married+tr_IslamNegative*married+
                 dm1+dm2+dm3_edu+unempl+married
               +april+may+june+july, 
               data=data, family = binomial(link = "probit"))
coeftest(model4.16)


#Unemployed
model4.9<-polr(dv_UnderstandViolence~tr_IslamPositive*unempl+
                 tr_IslamNegative*unempl+dm3_edu+unempl+married+dm1+dm2
               +april+may+june+july, 
               data=data, method=c('probit'))
coeftest(model4.9)

model4.10<-polr(dv_Muslimidentity~tr_IslamPositive*unempl+tr_IslamNegative*unempl+
                 dm3_edu+unempl+married+dm1+dm2
               +april+may+june+july, 
               data=data, method=c('probit'))
coeftest(model4.10)

model4.17<-glm(dv_USfavorable~tr_IslamPositive*unempl+tr_IslamNegative*unempl+
                 dm1+dm2+dm3_edu+unempl+married
               +april+may+june+july, 
               data=data, family = binomial(link = "probit"))
coeftest(model4.17)

#Educated
model4.11<-polr(dv_UnderstandViolence~tr_IslamPositive*dm3_edu+
                 tr_IslamNegative*dm3_edu+dm3_edu+unempl+married+dm1+dm2
               +april+may+june+july, 
               data=data, method=c('probit'))
coeftest(model4.11)

model4.12<-polr(dv_Muslimidentity~tr_IslamPositive*dm3_edu+tr_IslamNegative*dm3_edu+
                  dm3_edu+unempl+married+dm1+dm2
                +april+may+june+july, 
                data=data, method=c('probit'))
coeftest(model4.12)

model4.18<-glm(dv_USfavorable~tr_IslamPositive*dm3_edu+tr_IslamNegative*dm3_edu+
                 dm1+dm2+dm3_edu+unempl+married
               +april+may+june+july, 
               data=data, family = binomial(link = "probit"))
coeftest(model4.18)


#SI.6.1
stargazer(model4.2, model4.4, 
          model4.8, model4.10, model4.12, 
          star.cutoffs = c(0.05, 0.01, 0.001))

#SI.6.2
stargazer(model4.13, model4.14,
          model4.16, model4.17, model4.18, 
          star.cutoffs = c(0.05, 0.01, 0.001))

#SI.6.3
stargazer(model4.1, model4.3, 
          model4.7,model4.9, model4.11, 
          star.cutoffs = c(0.05, 0.01, 0.001))



#SI.6.4: Region
#1 is FBiH, 2 RS, 3 Brcko (excluded)
data$entity_dummy<-ifelse(data$entity==1, 1, NA)
data$entity_dummy<-ifelse(data$entity==2, 0, data$entity)

model5<-polr(dv_UnderstandViolence~tr_IslamPositive+tr_IslamNegative+
                april+may+june+july+dm1+dm2+dm3_edu+unempl+married, 
              data=data[data$entity_dummy==0,], method=c('probit'))
coeftest(model5)

model5.1<-polr(dv_UnderstandViolence~tr_IslamPositive+tr_IslamNegative+
                april+may+june+july+dm1+dm2+dm3_edu+unempl+married, 
              data=data[data$entity_dummy==1,], method=c('probit'))
coeftest(model5.1)


model5.2<-polr(dv_Muslimidentity~tr_IslamPositive+tr_IslamNegative+
               april+may+june+july+dm1+dm2+dm3_edu+unempl+married, 
             data=data[data$entity_dummy==0,], method=c('probit'))
coeftest(model5.2)

model5.3<-polr(dv_Muslimidentity~tr_IslamPositive+tr_IslamNegative+
                 april+may+june+july+dm1+dm2+dm3_edu+unempl+married, 
               data=data[data$entity_dummy==1,], method=c('probit'))
coeftest(model5.3)


model5.4<-glm(dv_USfavorable~tr_IslamPositive+tr_IslamNegative+
                april+may+june+july+dm1+dm2+dm3_edu+unempl+married, 
              data=data[data$entity_dummy==0,], family = binomial(link = "probit"))
coeftest(model5.4)

model5.5<-glm(dv_USfavorable~tr_IslamPositive+tr_IslamNegative+
                april+may+june+july+dm1+dm2+dm3_edu+unempl+married, 
              data=data[data$entity_dummy==1,], family = binomial(link = "probit"))
coeftest(model5.5)


stargazer(model5.2, model5.3, model5.4, model5.5, model5, model5.1,
          star.cutoffs = c(0.05, 0.01, 0.001))



